/*
 * RayTracer.cpp
 * This class contains all the functionality for performing recursive ray tracing of the scene
 * - stores all objects in the scene
 * - contains the Grid array which divides the objects polys
 * - supports triangles, planes, cylinders and spheres
 *
 * @version 	1.00 2 Feb 2010
 * @author 	Eric Stock
 */

#define BSP
#include "RayTracer.h"
#include "Material.h"

using namespace std;

int RayTracer::noSphereIntersections = 0;

const char* RayTracer::strSceneFile = "Scene/cornell.dat";

const flt RayTracer::smallValue = 0.001f;

/** Constructor RayTracer
 *	
 *	Initializes pointers to null
 *  Initializes standard type variables to 0
 *
 */
RayTracer::RayTracer()

: m_treeBSP(1024, 1024)
, m_pLights(NULL)
, m_pAreaLightBoxes(NULL)
, m_pSpheres(NULL)
, m_pPlanes(NULL)
, m_pImageData(NULL)
, m_pSamples(NULL)
, m_fCameraCircleRadius(1.0)
, m_fCameraCircleAngle(0.0)
, m_iScreenHeight(0)
, m_iScreenWidth(0)
, m_iNumberLights(0)
, m_iNumberSpheres(0)
, m_iNumberPlanes(0)
{}

RayTracer::~RayTracer() {
   
   delete [] m_pSamples;
   delete [] m_pLights;
   delete [] m_pSpheres;
   delete [] m_pPlanes;
   delete [] m_pImageData;
   delete [] m_arrPrimitives;
}

/** Method TraceRay
 *	
 *	This is the recursive call which is called to intersect 
 *		-primary rays from origin
 *		-reflected rays (often up to a limit of recursive reflective calls)
 *		-refracted rays
 *	Performs lighting calculations for all intersected objects
 *		-ambient
 *		-diffuse
 *		-specular
 *
 *	@param int count - counter for number of recursively generated reflective rays
 *	@param Vector* origin - ray origin
 *	@param Vector* direct - ray direction
 *	@param Colour& colour - the colour of the current pixel being calculated through ray tracing algorithm
 *	@param flt colourAmount - coefficient of how much this ray can contribute to the pixel's colour (0.0 to 1.0)
 *	@param flt n1 - refractive index for material ray is travelling through
 *
 */
void RayTracer::TraceRay(int p_iCount, math::Vertex& p_vOrigin, math::Vector& p_vDirection, Colour& p_Colour, flt p_fColourAmount, flt p_fN1) 
{	
	Vertex closestPt;	/** point of intersection for closest object */				
	flt u = -1.0;
	flt v = -1.0;
    Vector normal;
   
	flt fEnter = 0.0f;
	flt fExit = 0.0f;

	bool bHitLight = m_pAreaLightBoxes[0]->IntersectRay(p_vOrigin, p_vDirection, 0, 100000.0f, fEnter, fExit);

	// intersect the BSP tree which contains spheres
	CollisionInfo info;
	flt fMaxDistance = 1000000.0f;
	if(bHitLight)
	{
		fMaxDistance = fEnter;
	}
	
	info = m_treeBSP.RecursiveRayTreeTraversalB(p_vOrigin, p_vDirection, fMaxDistance);
	
	if(info.pObject)
	{
		Primitive* pPrimitive = info.pObject;
		/** since we have an intersected object, we must perform shading for each light in the scene */
		for(int m=0; m < m_iNumberLights; m++) {	
			
			if(info.pObject->GetPrimitiveType() == Primitive::AP_Sphere) {
				
				flt refract = pPrimitive->m_Material.m_fRefract;

				Vector surface = pPrimitive->GetNormalAtIntersection(info.vIntersectPt);
				surface *= (flt)info.insideSphere;
				
				if(( pPrimitive->m_Material.m_fRefract > 0.0f || pPrimitive->m_Material.m_fReflect > 0.0f) && m == 0 && p_iCount < 4) {
					p_iCount++;
					flt dotProd = (surface[0] * p_vDirection[0]) + (surface[1] * p_vDirection[1]) + (surface[2] * p_vDirection[2]);	
					
					if(pPrimitive->m_Material.m_fReflect > 0.0f && p_fN1 < 1.32f) {
						Vector reflect(p_vDirection[0] - (2.0f * dotProd * surface[0]), p_vDirection[1] - (2.0f * dotProd * surface[1]), p_vDirection[2] - (2.0f * dotProd * surface[2]));
						reflect.Normalize();
						TraceRay(p_iCount, info.vIntersectPt, reflect, p_Colour, p_fColourAmount * pPrimitive->m_Material.m_fReflect, p_fN1);
					} 
					
					if(pPrimitive->m_Material.m_fRefract > 0.0f)
					{	
						flt n2 = 1.0f;
	            
						if(info.insideSphere ==  1)
						{
							n2 = pPrimitive->m_Material.m_fRefractIndex;			/** the sphere is full of water */
						}
						
						flt sin2Trans = pow(p_fN1/n2, 2.0f) * (1.0f - dotProd * dotProd);
						if(sin2Trans < 1.0f)
						{	
							math::Vector vTrans(p_fN1/n2 * p_vDirection[0] - (p_fN1/n2 * dotProd + sqrt(1.0f-sin2Trans)) * surface[0],
												p_fN1/n2 * p_vDirection[1] - (p_fN1/n2 * dotProd + sqrt(1.0f-sin2Trans)) * surface[1],
												p_fN1/n2 * p_vDirection[2] - (p_fN1/n2 * dotProd + sqrt(1.0f-sin2Trans)) * surface[2]);
							vTrans.Normalize();
							TraceRay(p_iCount, info.vIntersectPt, vTrans, p_Colour, p_fColourAmount * pPrimitive->m_Material.m_fRefract, n2);
						}
					}  
				}
				flt fLightCoeff = m_pLights[m]->GetLightContribution(info.vIntersectPt, &m_treeBSP);
				if(fLightCoeff > 0)
				{
					p_Colour += pPrimitive->Shade(*(m_pLights[m]), info.vIntersectPt,  p_vDirection, surface) * fLightCoeff;
				}	
				p_Colour += pPrimitive->GetAmbientColour(m_ColourAmbient, p_fColourAmount);
			}
			else if(info.pObject->GetPrimitiveType() == Primitive::AP_Plane) 
			{	
				Vector& vNormal = pPrimitive->GetNormalAtIntersection(info.vIntersectPt);

				if(pPrimitive->m_Material.m_fReflect > 0.0f && m==0 && p_iCount < 4)
				{
					Vector reflect(p_vDirection[0]-(2.0f * info.planeDotProduct * vNormal.data[0]),
							p_vDirection[1]-(2.0f * info.planeDotProduct * vNormal.data[1]),
							p_vDirection[2]-(2.0f * info.planeDotProduct * vNormal.data[2]));
					reflect.Normalize();
					
					TraceRay(++p_iCount, info.vIntersectPt, reflect, p_Colour, p_fColourAmount * pPrimitive->m_Material.m_fReflect, p_fN1);
				}
				
				flt fLightCoeff = m_pLights[m]->GetLightContribution(info.vIntersectPt, &m_treeBSP);
				if(fLightCoeff > 0)
				{
					p_Colour += pPrimitive->Shade(*(m_pLights[m]), info.vIntersectPt, p_vDirection, vNormal) * fLightCoeff;
				}
				p_Colour += pPrimitive->GetAmbientColour(m_ColourAmbient, p_fColourAmount);
			}
		}
	}	
	else if(bHitLight)
	{
		p_Colour.R += 1.0f * p_fColourAmount;
		p_Colour.G += 1.0f * p_fColourAmount;
		p_Colour.B += 1.0f * p_fColourAmount;
	}
}


void  RayTracer::LaunchThread(void* p_pParam) {

      ThreadArgs* pArgs = (ThreadArgs*)p_pParam;
	  RayTracer * pTracer = pArgs->tPtr;

      pTracer->IntersectScene(pArgs->scanLineId);       // now call the true entry-point-function

      // A thread terminates automatically if it completes execution,
      // or it can terminate itself with a call to _endthread().

      //return 1;          // the thread exit code

}

void RayTracer::NextFrame()				//repositions items to display the next frame
{
	m_pSpheres[1].m_vCenter.data[0] += 0.2;
	m_pSpheres[0].m_vCenter.data[1] -= 1.0;
	m_pSpheres[3].m_vCenter.data[2] += 2.0;
	m_pSpheres[4].m_vCenter.data[0] += 0.3;
	m_pSpheres[4].m_vCenter.data[0] -= 0.3;

	m_pSpheres[5].m_vCenter.data[0] += 0.2;
	m_pSpheres[6].m_vCenter.data[1] -= 0.3;
	m_pSpheres[7].m_vCenter.data[2] += 2.0;
}


void RayTracer::PhotonMapScene()
{
	m_pPhotons = new Photon[2000];		//store a max of 1000 photons
	int iTotalEmittedPhotons = 0;
	
	for(int i = 0; i < m_iNumberLights; i++)
	{
		flt x,y,z;
		int iEmittedPhotons = 0;
		Colour colour(m_pLights[i]->R, m_pLights[i]->G, m_pLights[i]->B);
		while(iEmittedPhotons < 1000)
		{
			do 
			{
				x = 2.0f * double((rand())/RAND_MAX) - 1.0f;
				y = 2.0f * double((rand())/RAND_MAX) - 1.0f;
				z = 2.0f * double((rand())/RAND_MAX) - 1.0f;
			} while(x*x + y*y + z*z > 1);
			
			math::Vector vDirection(x, y, z);
			vDirection.Normalize();
			
			//TraceRay(..pass in photons
			iEmittedPhotons++;
			iTotalEmittedPhotons++;
		}
	}
}

void RayTracer::TestRandomHemisphereVector()
{
	for(int i = 0; i < 100 ; i++)
	{
		math::Vector vTest = GetRandomNormalInHemisphere(math::Vector(0.0f, 0.0f, 0.5f));
		printf("x is %f y is %f z is %f\n",vTest[0],vTest[1],vTest[2]);
	}
}

/** Method IntersectScene
 *	
 *	For every pixel in the virtual screen, fire the primary rays through the pixel
 *  Note, currently we are using 4X anti-aliasing
 *  
 */
void RayTracer::IntersectScene(int p_iScanLineId)      {

	int index;
	Colour colour;
    
	flt refractiveIndex = 1.0;       /** refractive index for air */
    for(int i=0; i < m_iScreenHeight; i++) {
      //if(i%3 != p_iScanLineId) continue;
	  //if(iScanLineId > 0) continue;
	  //printf("row is %d\n",i);	
	  for(int j=0 ;j < m_iScreenWidth; j++) {
			int colourIndex = j * 3 + (i * m_iScreenWidth * 3);
            index = j + (i * m_iScreenWidth);
			colour.R = 0.0;
			colour.G = 0.0;
			colour.B = 0.0;
            
			for(int m=0;m<1;m++) {
				TraceRay(0, m_vOrigin, m_pSamples[index].rays[m], colour, 1.0, refractiveIndex);
			}
			//imageData[index].x = j;
			//imageData[index].y = i;
			//imageData[index].colour[0]=colour.R/4.0f;
			//imageData[index].colour[1]=colour.G/4.0f;
			//imageData[index].colour[2]=colour.B/4.0f;	
			
			m_pImageData[colourIndex]=colour.R;
			m_pImageData[colourIndex+1]=colour.G;
			m_pImageData[colourIndex+2]=colour.B;
		}
	}
	//printf("thread %d finishing \n",iScanLineId);
}

/** Method IntitializeScene
 *	
 *	arses the scene input file by opening the file and reading the data into variables
 *  
 *  
 */	
void RayTracer::IntitializeScene() {
	int i;
	int j;
	
	FILE * in_file;
	j=1;
	/** open the scene file */
	if((in_file = fopen(strSceneFile,"r"))==NULL) {	
			printf("file open failed\n");
	} else {
		/** parse the input file */
		fscanf(in_file,"%d%d", &m_iScreenWidth, &m_iScreenHeight);			/** width , height */
		ReadComment(in_file);								/** '/' char seperates information in file*/
		fscanf(in_file,"%f%f%f", &m_vOrigin.data[0], &m_vOrigin.data[1], &m_vOrigin.data[2]);		/** camera origin */
		ReadComment(in_file);
		math::Vertex vBoxMin, vBoxMax;

		fscanf(in_file,"%f%f%f", &vBoxMin.data[0], &vBoxMin.data[1], &vBoxMin.data[2]);			//starting min and max extents for bsp tree
		ReadComment(in_file);
		fscanf(in_file,"%f%f%f", &vBoxMax.data[0], &vBoxMax.data[1], &vBoxMax.data[2]);
		ReadComment(in_file);
		fscanf(in_file,"%f%f%f", &m_ColourAmbient.R, &m_ColourAmbient.G, &m_ColourAmbient.B);	/** ambient light */
		ReadComment(in_file);
		fscanf(in_file,"%d", &m_iNumberLights);							/** number of lights */
		m_pLights = new Light*[m_iNumberLights];
		m_pAreaLightBoxes = new Tree::Box*[m_iNumberLights];	
		ReadComment(in_file);
		
		/** initialize each light */
		for(i = 0; i < m_iNumberLights; i++) {						
			AreaLight* pLight = new AreaLight();
			m_pAreaLightBoxes[i] = new Tree::Box();
			m_pLights[i] = pLight;
			pLight->m_fLength = 10;
			
			//m_pLights[i] = new AreaLight();
			fscanf(in_file,"%f%f%f", &m_pLights[i]->m_vPosition[0], &m_pLights[i]->m_vPosition[1], &m_pLights[i]->m_vPosition[2]);
			ReadComment(in_file);
			fscanf(in_file,"%f%f%f", &m_pLights[i]->R, &m_pLights[i]->G, &m_pLights[i]->B);
			ReadComment(in_file);
			m_pAreaLightBoxes[i]->m_vMin.data[0] = m_pLights[i]->m_vPosition[0] - pLight->m_fLength/2.0f;
			m_pAreaLightBoxes[i]->m_vMin.data[1] = m_pLights[i]->m_vPosition[1];// - 0.1f;
			m_pAreaLightBoxes[i]->m_vMin.data[2] = m_pLights[i]-> m_vPosition[2] - pLight->m_fLength/2.0f;

			m_pAreaLightBoxes[i]->m_vMax.data[0] = m_pLights[i]->m_vPosition[0] + pLight->m_fLength/2.0f;
			m_pAreaLightBoxes[i]->m_vMax.data[1] = m_pLights[i]->m_vPosition[1];// + 0.1f;
			m_pAreaLightBoxes[i]->m_vMax.data[2] = m_pLights[i]->m_vPosition[2] + pLight->m_fLength/2.0f;
		}
		fscanf(in_file,"%d",&m_iNumberSpheres);		
		ReadComment(in_file);
		if(m_iNumberSpheres > 0)
		{
			m_pSpheres = new Sphere2[m_iNumberSpheres];
			for(i = 0; i < m_iNumberSpheres; i++) {
				fscanf(in_file,"%f%f%f%f", &m_pSpheres[i].m_vCenter[0], &m_pSpheres[i].m_vCenter[1], &m_pSpheres[i].m_vCenter[2], &m_pSpheres[i].m_fRadius);
				m_pSpheres[i].m_fRadiusSquared = m_pSpheres[i].m_fRadius * m_pSpheres[i].m_fRadius;
				ReadComment(in_file);
				fscanf(in_file,"%f%f%f", &m_pSpheres[i].m_Material.m_Colour.R, &m_pSpheres[i].m_Material.m_Colour.G, &m_pSpheres[i].m_Material.m_Colour.B);
				ReadComment(in_file);
				fscanf(in_file,"%f%f%f", &m_pSpheres[i].m_Material.m_fDiffuse, &m_pSpheres[i].m_Material.m_fReflect, &m_pSpheres[i].m_Material.m_fRefract);
				ReadComment(in_file);
				fscanf(in_file,"%f%f%f", &m_pSpheres[i].m_Material.m_fSpecular, &m_pSpheres[i].m_Material.m_fSpecularExp, &m_pSpheres[i].m_Material.m_fRefractIndex);
				ReadComment(in_file);
			} 
		}

		fscanf(in_file, "%d", &m_iNumberPlanes);			
		ReadComment(in_file);
		m_pPlanes = new Plane2[m_iNumberPlanes];
 		for(i = 0; i < m_iNumberPlanes; i++) {
			fscanf(in_file,"%f%f%f%f",&m_pPlanes[i].m_vNormal.data[0], &m_pPlanes[i].m_vNormal.data[1], &m_pPlanes[i].m_vNormal.data[2], &m_pPlanes[i].m_fD);
			m_pPlanes[i].m_vNormal.Normalize();
			ReadComment(in_file);
			fscanf(in_file,"%f%f%f", &m_pPlanes[i].m_Material.m_Colour.R, &m_pPlanes[i].m_Material.m_Colour.G, &m_pPlanes[i].m_Material.m_Colour.B);
			ReadComment(in_file);
			fscanf(in_file,"%f%f%f", &m_pPlanes[i].m_Material.m_fDiffuse, &m_pPlanes[i].m_Material.m_fReflect, &m_pPlanes[i].m_Material.m_fRefract);
			ReadComment(in_file);
			fscanf(in_file,"%f%f", &m_pPlanes[i].m_Material.m_fSpecular, &m_pPlanes[i].m_Material.m_fSpecularExp);
			ReadComment(in_file);
		}

		fclose(in_file);
	
		int iTotalPrimitives = m_iNumberSpheres + m_iNumberPlanes;
		
		m_arrPrimitives = new Primitive*[iTotalPrimitives];
		
		for(int i = 0; i < m_iNumberSpheres; i++)
		{
			m_arrPrimitives[i] = (Primitive*)&m_pSpheres[i];
		}

		for(int i = 0; i < m_iNumberPlanes; i++)
		{
			m_arrPrimitives[m_iNumberSpheres + i] = (Primitive*)&m_pPlanes[i];
		}
		
		m_treeBSP.BuildTree(m_arrPrimitives , iTotalPrimitives, m_pAreaLightBoxes, m_iNumberLights, vBoxMin, vBoxMax);
		
		/** these are used for rotating the camera during animation */
		m_fCameraCircleRadius = m_vOrigin.data[2];
		m_fCameraCircleAngle = 4.71f;
		
	    printf("Scene sucessfully read in \n");
	}
}

/** Method UpdateCameraPos
 *	
 *	Divides a complete rotation around the scene to 300 frames
 *  -frame 0 is the camera's position in the .dat file
 *
 *  @param int frame - the frame number from 0 to 299
 *  
 */	
void RayTracer::UpdateCameraPos(int p_iFrame) {
     
	flt tempAngle = m_fCameraCircleAngle + p_iFrame * (2.0f * 3.14f)/300.0f;
    m_vOrigin.data[0] = m_fCameraCircleRadius * cos(tempAngle);
    m_vOrigin.data[2] = m_fCameraCircleRadius + m_fCameraCircleRadius * sin(tempAngle);   
}

		
/** Method ReadComment
 *	
 *	Parses the delimiter '/' in the input file
 *
 *  @param FILE * in_file - pointer to file being read in
 *  
 */	
void RayTracer::ReadComment(FILE * p_pFile) {
	int j;
	char c;
	j=1;
	while(j) {
		fscanf(p_pFile,"%c",&c);			//read a character at a time
		if(c=='/')break;						//until '/'
	}
}

math::Vector RayTracer::GetRandomNormalInHemisphere(math::Vector& p_vNormal)
{
	flt fRand1 = double((double)(rand())/(double)RAND_MAX);
	flt fRand2 = double((double)(rand())/(double)RAND_MAX);

	flt fPolarAngle = acos(sqrt(fRand1));
	flt fAzimuthalAngle = 2.0f * M_PI * fRand2;

	//convert the spherical coords to cartesian coordinates
	flt x = sin(fPolarAngle) * cos(fAzimuthalAngle);
	flt z = sin(fPolarAngle)  * sin(fAzimuthalAngle);
	flt y = cos(fPolarAngle);
	
	//create 2 vectors to form a basis with the normal
	math::Vector vBasis1;
	
	//create 1st vector by crossing with standard basis vector with least magniture
	
	if(p_vNormal[0] <= p_vNormal[1] && p_vNormal[0] <= p_vNormal[2])
	{
		math::Vector vXaxis(1.0f, 0.0f, 0.0f);
		vBasis1 = math::Vector::Cross(p_vNormal, vXaxis);
	}
	else if(p_vNormal[1] <= p_vNormal[0] && p_vNormal[1] <= p_vNormal[2])
	{
		math::Vector vYaxis(0.0f, 1.0f, 0.0f);
		vBasis1 = math::Vector::Cross(p_vNormal, vYaxis);
	}
	else if(p_vNormal[2] <= p_vNormal[0] && p_vNormal[2] <= p_vNormal[1])
	{
		math::Vector vZaxis(0.0f, 0.0f, 1.0f);
		vBasis1 = math::Vector::Cross(p_vNormal, vZaxis);
	}

	//create 2nd basis vector by crossing basis vector 1 with parameter normal vector

	math::Vector vBasis2 = math::Vector::Cross(vBasis1, p_vNormal);
	
	p_vNormal = p_vNormal * y;
	vBasis1 = vBasis1 * x;
	vBasis2 = vBasis2 * z;

	math::Vector randVector(p_vNormal[0] + vBasis1[0] + vBasis2[0], 
							p_vNormal[1] + vBasis1[1] + vBasis2[1],
							p_vNormal[2] + vBasis1[2] + vBasis2[2]);
	
	return randVector;

	//rotate 90 in x axis

	
	/*
	
	//calc angle between world up vector and normal vector
	flt angle = acos((p_vNormal[1]));
	//angle *= -1.0f;
	math::Vector rotAxis = math::Vector::Cross(p_vNormal, math::Vector(0.0f, 1.0f, 0.0f));

	//create matrix for rotation about rotation matrix with calculated angle

	math::Matrix rotMatrix;
	
	flt t = 1.0f - cos(angle);
	flt c = cos(angle);
	flt s = sin(angle);

	rotMatrix[0] = t*rotAxis[0]*rotAxis[0] + c;
	rotMatrix[1] = t*rotAxis[0]*rotAxis[1] - s*rotAxis[2];
	rotMatrix[2] = t*rotAxis[0]*rotAxis[2] + s * rotAxis[1];
	rotMatrix[3] = 0;

	rotMatrix[4] = t*rotAxis[0]*rotAxis[1] + s*rotAxis[2];
	rotMatrix[5] = t*rotAxis[1]*rotAxis[1] + c;
	rotMatrix[6] = t*rotAxis[1]*rotAxis[2] - s * rotAxis[0];
	rotMatrix[7] = 0;

	rotMatrix[8] = t*rotAxis[0]*rotAxis[2] - s*rotAxis[1];
	rotMatrix[9] = t*rotAxis[1]*rotAxis[2] + s*rotAxis[0];
	rotMatrix[10] = t*rotAxis[2]*rotAxis[2] + c;
	rotMatrix[11] = 0;

	rotMatrix[12] = 0;
	rotMatrix[13] = 0;
	rotMatrix[14] = 0;
	rotMatrix[15] = 1;

	math::Vector randVector;

	randVector[0] = x * rotMatrix[0] + y * rotMatrix[4] + z * rotMatrix[8];
	randVector[1] = x * rotMatrix[1] + y * rotMatrix[5] + z * rotMatrix[9];
	randVector[2] = x * rotMatrix[2] + y * rotMatrix[6] + z * rotMatrix[10];

	return randVector;
	*/


}

/** Method GenScreenDirections
 *	
 *	This function generates all the normalized primary direction rays
 *	Takes into consideration the camera field of view
 *  
 */	
void RayTracer::GenScreenDirections() {
			
	int numberPixels;
	int index;
	Vector temp;			
	Matrix trans;
    numberPixels = m_iScreenWidth * m_iScreenHeight;		
 	m_pSamples = new Samples[numberPixels];
  	m_pImageData = new Pics[numberPixels * 3];		
	
	Vector targetPosition;

	/** here is our hard coded target for the camera 
	*	to-do add the ability to set this in the .dat file
	*/
	targetPosition[0] = 0.0;
	targetPosition[1] = 0.0;
	targetPosition[2] = 0.0;
	
	/** here the vector space for the camera is built using the look at as the z axis
	*	the cross between the z and world up as the x axis and the cross between 
	*/

	Vector zAxis;									 /** zAxis vector represents the look at */
    zAxis[0] = targetPosition[0] - m_vOrigin[0];                
	zAxis[1] = targetPosition[1] - m_vOrigin[1];  
	zAxis[2] = targetPosition[2] - m_vOrigin[2];  
    zAxis.Normalize();
	
    Vector up;										/** world up vector */
	up[0] = 0.0;
	up[1] = 1.0;
	up[2] = 0.0;
	
	Vector xAxis = Vector::Cross(up, zAxis);            /** xAxis */
	xAxis.Normalize();
    Vector yAxis = Vector::Cross(zAxis, xAxis);		   /** yAxis */
	yAxis.Normalize();

    /** build the transform matrix - combined rotation and translation */
    trans.cell[0] = xAxis[0]; trans.cell[1] = xAxis[1]; trans.cell[2] = xAxis[2];
    trans.cell[4] = yAxis[0]; trans.cell[5] = yAxis[1]; trans.cell[6] = yAxis[2];
    trans.cell[8] = zAxis[0]; trans.cell[9] = zAxis[1]; trans.cell[10] = zAxis[2];
    
    //flt tanFov = 1.33;								/**  deg field of view */
	flt tanFov = 0.414f;								/** 45 deg field of view */
	flt leftScreen = -1.0 * tanFov;
	flt rightScreen = 1.0 * tanFov;
	flt bottomScreen = -1.0 * (flt)m_iScreenHeight / (flt)m_iScreenWidth *  tanFov;
    flt topScreen = 1.0 * (flt)m_iScreenHeight / (flt)m_iScreenWidth * tanFov; 
    
	/** create virtual screen along z==1 plane */
    Vertex topLeft(leftScreen, topScreen, 1.0);
    Vertex topRight(rightScreen, topScreen, 1.0);
    Vertex bottomLeft(leftScreen, bottomScreen, 1.0);
    Vertex bottomRight(rightScreen, bottomScreen, 1.0);
    
	/** transform this screen into world space */
	topLeft = trans.TransformVert(topLeft);
    topRight = trans.TransformVert(topRight);
    bottomLeft = trans.TransformVert(bottomLeft);
    bottomRight = trans.TransformVert(bottomRight);
    
    /** dX is the vector which must be added as the ray moves across the virtual screen */
	Vector dX;
    dX[0] = (topRight[0] - topLeft[0])/(flt)m_iScreenWidth;
	dX[1] = (topRight[1] - topLeft[1])/(flt)m_iScreenWidth; 
	dX[2] = (topRight[2] - topLeft[2])/(flt)m_iScreenWidth;
    
	/** dY is the vector which must be added as the ray moves up along the virtual screen */
	Vector dY;
    dY[0] = (topRight[0] - bottomRight[0]) / (flt)m_iScreenHeight; 
	dY[1] = (topRight[1] - bottomRight[1]) / (flt)m_iScreenHeight; 
	dY[2] = (topRight[2] - bottomRight[2]) / (flt)m_iScreenHeight;
    
	/** here are the 3 multi samples which go along with bottomLeft ray */
	Vertex bottomLeftMulti1(leftScreen + dX[0]/2.0f, bottomScreen, 1.0);                    /** sample slightly to left of bottom left */ 
    Vertex bottomLeftMulti2( leftScreen, bottomScreen + dY[1]/2.0, 1.0);                    /** sample slightly above of bottom left */ 
    Vertex bottomLeftMulti3(leftScreen + dX[0]/2.0, bottomScreen + dY[1]/2.0, 1.0);        /** sample slightly to right of bottom left */
 
    /** transform the 3 multi samples */
    bottomLeftMulti1 = trans.TransformVert(bottomLeftMulti1);
    bottomLeftMulti2 = trans.TransformVert(bottomLeftMulti2);
    bottomLeftMulti3 = trans.TransformVert(bottomLeftMulti3);
    
	/** loop through each pixel in the virtual screen and calculate the ray direction */
    for(int i=0; i < m_iScreenHeight; i++) {
        for(int j=0; j < m_iScreenWidth; j++) {
			
			index = (j + i * m_iScreenWidth);
            
			m_pSamples[index].rays[0][0] = bottomLeft[0] + j * dX[0];
			m_pSamples[index].rays[0][1] = bottomLeft[1] + i * dY[1];
			m_pSamples[index].rays[0][2] = bottomLeft[2];
			m_pSamples[index].rays[0].Normalize();

			m_pSamples[index].rays[1][0] = bottomLeftMulti1[0] + j * dX[0];
			m_pSamples[index].rays[1][1] = bottomLeftMulti1[1] + i * dY[1];
			m_pSamples[index].rays[1][2] = bottomLeftMulti1[2] + j * dX[2];
			m_pSamples[index].rays[1].Normalize(); 
			
			m_pSamples[index].rays[2][0] = bottomLeftMulti2[0] + j * dX[0];
			m_pSamples[index].rays[2][1] = bottomLeftMulti2[1] + i * dY[1];
			m_pSamples[index].rays[2][2] = bottomLeftMulti2[2] + j * dX[2];
			m_pSamples[index].rays[2].Normalize(); 

			m_pSamples[index].rays[3][0] = bottomLeftMulti3[0] + j * dX[0];
			m_pSamples[index].rays[3][1] = bottomLeftMulti3[1] + i * dY[1];
			m_pSamples[index].rays[3][2] = bottomLeftMulti3[2] + j * dX[2];
			m_pSamples[index].rays[3].Normalize(); 
		}
	}
}
	
